home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / spline_p.pro < prev    next >
Text File  |  1997-07-08  |  4KB  |  137 lines
  1. ; $Id: spline_p.pro,v 1.6 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1993-1997, Research Systems, Inc. All rights reserved.
  4. ;    Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;    SPLINE_P
  8. ;
  9. ; PURPOSE:
  10. ;    This procedure performs parameteric cubic spline interpolation.
  11. ;
  12. ; CATEGORY:
  13. ;    Interpolation - E1.
  14. ;
  15. ; CALLING SEQUENCE:
  16. ;    SPLINE_P, X, Y, Xr, Yr
  17. ;
  18. ; INPUTS:
  19. ;    X:      The abcissa vector (should be floating or double).
  20. ;    Y:      The vector of ordinate values corresponding to X.
  21. ;    Neither X or Y need be monotonic.
  22. ;
  23. ; KEYWORD PARAMETERS:
  24. ;    INTERVAL: The interval in XY space between interpolants. If
  25. ;          omitted, approximately 8 interpolants per XY segment
  26. ;          will result.
  27. ;    TAN0:      The tangent to the spline curve at X[0], Y[0]. If omitted,
  28. ;          the tangent is calculated to make the curvature of the
  29. ;          result zero at the beginning. This is a two element vector,
  30. ;          containing the X and Y components of the tangent.
  31. ;    TAN1:      The tangent to the spline curve at X[N-1], Y[N-1].If omitted,
  32. ;          the tangent is calculated to make the curvature of the
  33. ;          result zero at the end. This is a two element vector,
  34. ;          containing the X and Y components of the tangent.
  35. ;
  36. ; OUTPUTS:
  37. ;    XR:      The abcissa values of the interpolated function. This
  38. ;          may NOT be the same variable as either X or Y.
  39. ;    YR:      The ordinate values of the interpolated function. This
  40. ;          may NOT be the same variable as either X or Y.
  41. ;
  42. ; RESTRICTIONS:
  43. ;    X and Y should be floating or double.
  44. ;
  45. ; PROCEDURE:
  46. ;    Cubic spline interpolation with relaxed or clamped end conditions
  47. ;    as used in the Numerical Recipes.
  48. ;
  49. ;    This routine is both more general and faster than the
  50. ;    user's library function SPLINE. One call to SPLINE_P is equivalent
  51. ;    to two calls to SPLINE, as both the X and Y are interpolated with
  52. ;    splines. It is suited for interpolating between randomly 
  53. ;    placed points, and the abcissae    values need not be monotonic.
  54. ;    In addition, the end conditions may be optionally specified via
  55. ;    tangents.
  56. ;
  57. ; EXAMPLE:
  58. ;    The commands below show a typical use of SPLINE_P:
  59. ;      X = [0.,1,0,-1,0]      ;Abcissae for square with a vertical diagonal
  60. ;      Y = [0.,1,2,1,0]      ;Ordinates
  61. ;      SPLINE_P, X, Y, XR, YR  ;Interpolate with relaxed end conditions
  62. ;      PLOT, XR, YR          ;Show it
  63. ;
  64. ;     As above, but with setting both the beginning and end tangents:
  65. ;       SPLINE_P, X, Y, XR, YR, TAN0=[1,0], TAN1=[1,0]
  66. ;
  67. ;     This yields approximately 32 interpolants.
  68. ;
  69. ;     As above, but with setting the interval to 0.05, making more
  70. ;    interpolants, closer together:
  71. ;       SPLINE_P, X, Y, XR, YR, TAN0=[1,0], TAN1=[1,0], INTERVAL=0.05
  72. ;
  73. ;     This yields 116 interpolants and looks close to a circle.
  74. ;
  75. ; MODIFICATION HISTORY:
  76. ;    DMS, RSI.    August, 1993.    Written.
  77. ;    DMS, RSI.    Jan, 1994.  Modified to use NR_ spline routines.
  78. ;-
  79.  
  80. PRO SPLINE_P, x, y, xr, yr, $
  81.     INTERVAL=interval, TAN0=tan0, TAN1=tan1
  82.  
  83. n = n_elements(x)
  84. if n ne n_elements(y) then $
  85.     message,'X and Y must have the same number of points'
  86.  
  87. ni = n-1        ;Number of intervals
  88.  
  89. dx = x - shift(x,1)    ;Delta x and y
  90. dy = y - shift(y,1)    ;dx(i) = x(i) - x(i-1)
  91. dx[0] = 0.
  92. dy[0] = 0.
  93.  
  94. t = sqrt(dx^2 + dy^2)    ;interpoint Distance
  95. ni = n-1
  96. big = 2.0e30
  97.  
  98. ; Default interval = approx 8 points per interval....
  99. if n_elements(interval) le 0 then interval = total(t) / (8*ni)
  100.  
  101. r = ceil(t/interval)        ;# of elements in each interval
  102. nr = long(total(r))        ;# of elements in result
  103.  
  104. tt = fltarr(nr+1, /nozero)
  105. j = 0L
  106.  
  107. for int = 0, ni-1 do begin    ;Each interval
  108.     i1 = int+1
  109.     nn = r[i1]            ;# pnts in this interval
  110.     tt[j] = t[i1] / nn * findgen(nn) + t[int]
  111.     t[i1] = t[i1] + t[int]
  112.     j = j + nn
  113.     endfor
  114. tt[nr] = t[int]
  115.  
  116. ; Use end tangents, or use relaxed condition.
  117. if n_elements(tan0) ge 2 then begin    ;Clamped on left?
  118.     xp0 = tan0[0]
  119.     yp0 = tan0[1]
  120. endif else begin            ;Relaxed
  121.     xp0 = big
  122.     yp0 = big
  123. endelse
  124.  
  125. if n_elements(tan1) ge 2 then begin    ;Clamped on right?
  126.     xpn = tan1[0]
  127.     ypn = tan1[1]
  128. endif else begin
  129.     xpn = big
  130.     ypn = big
  131. endelse
  132.  
  133. ; Compute result & quit
  134. xr = SPL_INTERP(t,x, SPL_INIT(t,x, yp1 = xp0, ypn = xpn), tt)
  135. yr = SPL_INTERP(t,y, SPL_INIT(t,y, yp1 = yp0, ypn = ypn), tt)
  136. end
  137.